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Despite recent advances in the speed of digital computers and in 
numerical algorithms for the solution of differential equations, the 
evaluation of the dynamic response of nonlinear circuits is still too 
slow to permit Monte Carlo tolerance analysis. However, the per- 
formance of many nonlinear circuits can be evaluated on a static 
basis. One example is a D/A converter built with devices much faster 
than the converter's cycle time. Algorithms now exist that produce 
the static or equilibrium solution, of such networks in seconds. This 
paper deals with these algorithms and the associated techniques that 
have been embodied in a program for the Monte Carlo tolerance 
analysis of nonlinear, "dc," circuits. 

I. INTRODUCTION 

To date, most tolerance analysis of circuits has been in the fre- 
quency domain, as this series of articles indicates. The need for non- 
linear analysis arises not only for large signal circuits but also for 
small signal ac circuits where device model parameters vary with 
bias. Recent advances in the speed of digital computers and numerical 
algorithms have made possible the analysis of circuits with nonlinear 
behavior. Large signal, or time domain analysis of nonlinear circuits, 
however, is still such a comparatively slow process that Monte Carlo 
methods are out of the question. Enough algorithmic innovations have 
been achieved in the static analysis of nonlinear circuits that Monte 
Carlo methods can be applied to a wide class of nonlinear problems. 
DC in the sense used in this paper implies that the dynamic behavior 
of the nonlinear devices is fast in relation to the response of the rest 
of the circuit. There are, in general, three types of circuits that fall 
into this class. 

(i) Circuits that are essentially dc, such as operational amplifiers, 
power supplies, and those circuits used as examples in this paper. 
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(ii) Circuits that can be subdivided into dc and ac blocks where 
the nonlinear behavior of the circuit is present essentially only 
in the dc part. An example of this type of circuit is the Touch- 
Tone oscillator, the analysis of which is described in this 
series. 1 This class of circuits is very similar to the first. 
(in) Circuits designed for small signal applications may be analyzed 
in the frequency domain. Realistic modeling of devices, however, 
introduces changes in the small signal model parameters for 
different bias points. If the bias circuits vary randomly, then 
nonlinear dc analysis is required in the ac tolerance analysis 
loop. 

This paper deals with some of the methods required for the ef- 
ficient analysis of nonlinear circuits in a dc sense. The techniques 
and algorithms to be described have been embodied in a computer 
program which performs Monte Carlo tolerance analysis of nonlinear 
dc circuits as well as dc and transient analysis. Some of the more 
critical implementation aspects are described. 

II. DESIGN OF A NONLINEAR TOLERANCE ANALYSIS PROGRAM 

2.1 Problems Involved 

Until recently tolerance analysis, even nonlinear tolerance analysis, 
has been simple in concept. Circuits were manufactured using discrete 
components which generally had independent statistical behavior; 
transistors and diodes were expensive items and were used sparingly. 
Hence, a large amount of analysis could be done without a complex 
software package. 

Integrated circuits have opened a whole host of new problems in 
this area. Circuits designed today typically employ large numbers 
of transistors (since transistors are as cheap as resistors), posing many 
problems in their modeling, simulation, and solution. Probably the 
biggest dilemma in the design of a nonlinear tolerance analysis pro- 
gram is the minuscule past experience to draw upon as to what anal- 
ysis is required, what to do with the analysis results once they are 
obtained, and how to interpret them. 

Some of the special problems that arise in nonlinear tolerance anal- 
ysis are: 

(i) Parameters tend to be statistically correlated. This implies 
that many new output features have to be present in the soft- 
ware. 
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(it) There is a very wide range of circuit problems that will have 
to be solved. This poses special difficulties and restrictions on 
the methods of analysis, as will be seen in the following sections. 
(Hi) Statistical data for circuits manufactured today are not yet 
available or are available in limited quantities; 2 manufacturing 
processes vary from day to day and parameter correlations, 
aging data, etc., are all important to the analysis. 

2.2 Criteria To Be Met 

One of the most important criteria to be met in the design of a 
nonlinear tolerance analysis program is that it be easy to use. Some 
of the characteristics implied by this, both for the users of the pro- 
gram and for the writers of it, are: 

(i) The program must be humanly engineered to have a simple, 
clear, and easily learned input language. This applies not only 
to the network description but also to the description of sta- 
tistical data and to the command structure. The trend has been 
for engineers to personally use available computer tools rather 
than work through intermediaries, and the program itself 
should present as few obstacles as possible. In addition, the 
output capabilities must include data reduction schemes so 
that insight is gained at a glance. 
(ii) The program must be designed to be flexible enough so that it 
can be changed easily. Past experience has shown that a general- 
purpose analysis program will undergo many changes and, 
in fact, wall probably never reach a static condition. This means 
the program must be written in modular form, a characteristic 
that very often degrades efficiency. Modularity, however, 
implies ease of maintenance, upgrading, and debugging. 
(Hi) The program must be portable because of wide demand. This 
implies that it be written in some high level language, such as 
Fortran, with possibly a very small number of critical routines 
written in Assembly language for efficiency. 

In addition to ease of use, an important criterion is, of course, 
economy and reliability. The program must be very efficient if it is 
to be useful. Solution times for each statistical design must be mea- 
sured in seconds to make the analysis practical at all, and possibly 
in milliseconds if sophisticated features such as performance contours 
and large scale sensitivities are to be included. 3 Until recently, one 
would have been happy to get one solution to a nonlinear circuit; today 
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we are faced with obtaining hundreds and perhaps thousands of these 
solutions in a reasonable time. 

Various algorithms for solution of the nonlinear equations that 
arise from circuit simulations have been described in the literature; 
however, they are all basically variations on the Newton-Raphson 
scheme and are in general not suitable as they stand. The majority 
of the schemes converge in the order of tens of iterations, if they con- 
verge at all. Since solution times for each iteration are generally 
proportional to the cube of the number of variables, many iterations 
for each statistical solution preclude their use in a program such as 
this. Recent breakthroughs, however, in the Newton-Raphson solu- 
tion of nonlinear circuits and careful implementation of these methods 
permit solution times short enough (of the order of one second) so 
that meaningful tolerance analysis is possible. 

III. NUMERICAL ALGORITHMS AND THEIR IMPLEMENTATION 

3.1 Problem Formulation 

Given a network topology, there are various ways to write equa- 
tions describing the behavior of the network. Some examples are 
nodal equations, loop equations, Branin's 4 mixed mesh/cut-set equa- 
tions and the state-space formulation, which in the dc case has come 
to be known as the "normal form." 5 We have adopted the normal 
equation formulation for the following reasons : 

(i) The reduced set of equations produces, in general, a small 
system that has to be solved iteratively, by effectively separating 
the linear and nonlinear aspects of the problem. 
(ii) Implementation of various analysis techniques for equations 

in their normal form is straightforward. 
(Hi) The normal form of the equations can be generated extremely 
fast (see below) in a manner competitive with the most efficient 
sparse matrix techniques available today. 
Some other formulation (such as nodal equations), coupled with 
sparse matrix techniques may be more efficient for circuits containing 
a large number of junctions compared with the number of linear 
resistors. The formulation employed here, however, allows straight 
forward implementation of various convergence schemes such as the 
nonlinear transformation described in Section 3.3.3. In addition, it is 
felt that the normal formulation is the most practical for small-to- 
medium size circuits with a significant number of linear resistors 
(including those in the device models). The operational amplifier 
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example in this paper probably represents the limiting practical size 
of circuit for this formulation. Most of the methods described here are 
also applicable to any other equation formulation. 

3.2 Generation oj Equations in the Normal Form 

Consider any circuit consisting of current and voltage sources, 
diodes, transistors and resistors. From the network topology and net- 
work parameters, the large signal behavior of the network is charac- 
terized by a system of equations in the "normal form," viz., 

[Am + ram + [aw] = o (d 

where 

[VI is the vector of all device junction voltages and is the set of 
independent variables to be found; 

[U] is the vector of independent voltage and current sources in the 
network; 

[A] and [B] are coefficient matrices dependent on the network re- 
sistances; and 

[N(V)] is a vector of functions dependent on the nonlinear properties 
of the network. 
For computer simulation, the network devices are characterized by 
the Ebers-Moll model where the functional form of N(V) is 

N(V) = I D = I s [exp (6V) - 1] 

where Is is intercept current, and 6 is dependent on temperature. 
Notice that the formulation (1) isolates the linear part of the network 
from the nonlinear part and that the nonlinear behavior can be char- 
acterized by a vector quantity. 

The set of equations given by (1) must be generated for each sta- 
tistical design, since the [.A] and [B] matrices are dependent on 
resistor values. The implication of this, of course, is that a very fast 
algorithm is needed for equation generation. 

To completely characterize network behavior, another set of equa- 
tions is required which relates any requested network currents or 
voltages to the junction voltages found from (1). If [Z] is the vector 
of user-requested outputs (element currents and voltages), then 

[Z] = [C][V] + [D][U] (2) 

where again [C] and [D] are coefficient matrices dependent on resistor 
values and must be recalculated for each statistical design. 

The basis for the formulation of the A, B, C and D matrices is 



1178 



THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1971 



mainly topological in nature. From the network incidence matrix, [T], 
and the choice of a network tree as described below, the fundamental 
loop and cut-set matrices {[F] and [—F] T , respectively) are derived. 4 
The incidence matrix, T, has elements and ±1 and dimensions 
n X b for a network with (n + 1) nodes and b elements. Let I be a 
vector of element currents which is partitioned into tree branch currents, 
I„ and link currents, Ii . Then, 



[r][i] = [-*,*•] 



= o, 



I ( = FI, and V, = -F T V t . 



(3) 

(4) 



Equations (3) and (4) are expressions of Kirchoff's Current and Voltage 
Laws where d is the identity matrix and the vector of element voltages, 
V, is partitioned like I into V< and V! . [F] has dimensions n X (b — n), 
and its elements are and ±1. 

It is desired to place all voltage sources and device junctions in the 
tree and all current sources in links. Assuming this is possible, the 
columns of [T] are arranged in the following preference order (given 
also is the notation to be used for tree and link current and voltage) : 

/ V 



voltage sources 


Ib 


E 


device junctions 


h 


v D 


resistors — tree 


In 


V B 


—link 


la 


V a 


current sources 


J 


Vj 



The incidence matrix, [T], can be quickly reduced to the form 
[—/, F] by gaussian elimination, from which the fundamental loop 
matrix is denned. The reduction process favors the above top elements 
for inclusion in the tree. 

Now partition [F] as follows: 





G 


J 


E 


Fit 


F 12 


D 


F 2i 


F 22 


R 


F 3l 


F 3 2 



and from equation (4) write 

Id = F 21 I G 4- F 22 J, 



(5) 
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/r = F M U + F.J, (6) 

V = -F& - FZV D - FlV H . (7) 

Let [R] and \G] be defined as diagonal matrices whose elements are 
tree resistor values and link conductance values, respectively. Then, 
from equations (6) and (7), 



[~I* 




r iT 





1 


rv*~ 




v c _ 




L 


(T'J 


_i c _ 










= 






■F 2 \_ 


[V c ] + 




--*7. 



F 32 




E 

JJ 



+ 



F n 



_-fL o__i c _ 



v* 



or 



R- 1 

LSI 



-F n 
CT 1 



V* 

-I J 





L-F.TJ 



[V B ] + 







^3: 





E 

LJ. 



Define 



Then 



iT 1 -F„ 

LP. r i <T> . 






= [II]. 



= [H] 







Wd] + [H] 



L-Ff, 



F 3: 




E 

LJ. 



(8) 



(9) 



(10) 



Substituting I G into equation (5) results in equation (1), namely 

Id = N(V D ) = [A][V D ] + [B]\U] 
or 

AV + B\J - N(V) = 0. 

Also, from F and — F T , any network current or voltage can be extracted 
and the [Z] vector of equation (2) calculated. In general, 



[Z] = [C"][V] + [D'][U] + [P] 



Ir, 



(ID 



where C", D' and P are matrices whose elements are 0, ±1 obtained 
from selective rows and columns of F and — F T . Substituting equation 
(10) into equation (11) yields equation (2). 
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The time-consuming task in the calculation of the A, B, C, D matrices 
is in finding \ R as expressed by equation (10). This calculation 

requires one matrix inversion and two matrix multiplications. 

The number of rows and columns of H' 1 , as seen by equation (9), 
is equal to the number of resistors in the network. For a completely 
characterized network, including parasitics, this number can easily 
approach 80 or 100. The number of multiplications required to invert 
an nth order matrix is n 3 , or 10 6 for our example! Clearly, the imple- 
mentation of generating the equations of (1) and (2) is critical. 

Fortunately, the matrices that evolve from this formulation have 
special properties which can be advantageous. These are: 

(i) The fundamental loop matrix, F, contains only 0, ±1 entries 
and is sparse. From experience with a wide range of problems, 
this matrix is 20 to 50 percent dense (ratio of nonzero entries 
to total number of entries). Because of these properties, any 
matrix multiplication involving F or its partitions can be 
performed by additions and subtractions rather than by mul- 
tiplications. The operation of addition is at least 3 to 4 times 
faster than multiplication on most digital computers. Also, 
the number of additions and subtractions to be performed in 
multiplying F by A, where A is an Wi X n 2 matrix, is equal to 
the product of n 2 and the number of nonzero entries of the F 
matrix involved. 
(ii) In addition to F being sparse, other matrices as in equation (10), 
are sparse with predictably placed submatrices being identi- 
cally 0. 
(Hi) The special form of H' 1 of equation (9) allows profitable ap- 
plication of block inversion methods. This warrants further 
discussion. Consider equation (9) : 



IT 1 -F 31 



Hu Hi 
_H 21 H 2 



Both R' 1 and G' 1 are diagonal matrices and the dimension of 
H is n X n — (n R + n G ) 2 where n R is the number of tree re- 
sistors in the network and n G is the number of link resistances. 
A method of block inversion is chosen based on the min (n R , n G ). 

Case 1: n R ^ n Case 2: n Q > n R 

ffn = [R- 1 + F^GFlV, H 22 = [G' 1 + FjBF«r l , 
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II V2 = H u F ai O, II 21 = -H 22 FlR t 

H 21 = —GF-mHu, H l2 = RF 3 iH 22 , 

H 22 = G - GF^H^. H U =R + RF 3l H 2l . 

In case 1, the inversion of an n R X n R matrix (H n ) is required; for 
case 2 the matrix to be inverted is of order n G (#22). To invert H could 
take n 3 multiplications. The worst case for this formulation is n R = 
n G = n/2, requiring n 3 /8 multiplications. The calculation of the re- 
maining submatrices of H involve additions, subtractions or, at worst, 
multiplication by diagonal matrices. 

3.3 Solution of the dc Equations 

3.3.1 The Neivton-Raphson Method and Its Limitations 

The standard Newton-Raphson solution of equation (1) is obtained 
as follows: For an arbitrary estimate of V, say V*, equation (1) is not 
satisfied exactly, and there results: 

[A)[V k ] + [B]\U] + [N<y k )] = [R k ] (12) 

where the vector [R*] is termed the residual vector and is, in this case, 
a measure of the current imbalances in the circuit that result from 
insisting that [V] = [V*]. The superscript k refers to the iteration 
number. Successive estimates of [V] are formed as: 

[V* +I ] = [V] + P*] (13) 

where the step vector [P*] is obtained by solution of 

-[J(V*)1[P*] = [R*]. (14) 

In equation (14), [«7(V*)] is the Jacobian of the system (1), viz., 

[J(V k )] = [A] + [N'(V k )}. (15) 

The iterations are terminated whenever either the step vector [P] 
and/or the residual vector [R] is sufficiently close to zero. 

Straightforward application of the above iterative method to the 
system (1) results in several difficulties: 

(i) If the starting guess [V°] is not close to the solution, then there 
typically results exponential overflow in the nonlinear terms 
I s{e ' — 1), since the full Xewton-Raphson step may be too 
large in the positive direction. This may be overcome in several 
ways, the simplest probably being to reduce the step, 7 viz., 
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[V fc+1 ] = [V*] + SIP*] (16) 

where S is a step size, < S ^ 1 and is chosen so that 

22 (iO 2 < E (A*) 2 - (17) 

(ii) The method does not necessarily converge. Reliability can be 
greatly increased by the use of parameter-stepping techniques, 8 
as described in the next section. 

3.3.2 Increasing Reliability Through Source-Stepping 
Newton-Raphson methods typically converge whenever the starting 

guess is close to the solution point. This fact was first utilized by D. F. 
Davidenko 9 and subsequently by various authors. 10-12 Implementation 
of the method in nonlinear circuit analysis can be accomplished as 
follows 7 : For any system (1), one accurate solution is always known, 
namely, [V] = if [U] = 0. Hence, if the sources [U] are brought to 
their full value in small increments, convergence to each intermediate 
point is more likely. 

The strategy for stepping can take on various forms (see, for 
example, Ref. 12) ; the approach taken here involves source-stepping 
only when necessary. If no convergence is obtained after a fixed 
number of iterations with the sources on full, the sources are reduced 
to one-half their value and solution is again attempted. If necessary, 
the sources are progressively reduced until convergence is obtained at 
some intermediate source value, whereupon solution is again attempted 
with the sources on full. If convergence is not obtained, the sources 
are again reduced to a value midway between full and the last point 
at which convergence was obtained. The process is repeated until 
convergence at the full source value is obtained. 

3.3.3 Nonlinear Scalar Transjormation 

Reliability and speed of the Newton-Raphson method can also be 
greatly increased through the use of a nonlinear scalar transforma- 
tion of variables. The transformation is an extension 13 of the notion 
of "charge-state-variables" 14,15 based on a suitable definition of the 
"capacitance" of a junction, viz., 

q = f [1 + OKIs exp (6V)} dV (18) 

•In 

where (18) is a scalar equation. The potentials are considered to be 
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a function of q, i.e., V = V(q), and the Newton-Raphson iteration is 
performed in the g-space rather than the 7-space. 
This is accomplished by noting that 

[R(<?)1 = [R(7)J (19) 
and 

[J(q)) - [J(V)][S(V)r (20) 
where S(V) is a diagonal matrix with elements of the form, 

S it = 1 + BKJ a exp (OVi). (21) 

The Newton-Raphson step vector may be written in terms of q as 

P(«)] = -[^(3)] _, [R] (22) 

which results in 

Ptf )] = [S(V k )]\P(Y h )]. (23) 

What is required, then, is to transform the usual Newton-Raphson 
step vector, where now 

[q* +1 ] = [q*] + P(«*)]. (24) 

Since equation (23) represents a scalar transformation on the individual 
g, , equation (24) may be written in terms of the individual elements 
V { of the vector [V] as 

F* +1 + KJ Si exp (0F* +1 ) = Vl + KJ Si exp (07-) + 5,(7-)P,(7 •) 

(25) 

where 

S^VI) = 1 + 0KJ 8i exp (07*) 

and P,-(7*) is the ith element of the usual Newton-Raphson step 
in V. The transform parameter if, remains to be determined. 

Solution of equation (25) for each V** 1 gives the new estimate of 
the vector [V], given the standard Newton-Raphson step vector [P(7)], 
which may be obtained in the usual manner. Note that (25) is a scalar 
equation which may be solved very simply by Newton's method, as 
discussed further below. 

Empirical studies have shown that the speed of the iterative process 
is dependent on the value of K, which has the effect of "straightening 
out" the exponential characteristics at the expense of "warping" the 
linear parts of the solution space. It was determined that a good 
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choice for K is 

10 



Ki = 



10 2 < K < 10 4 



(26) 



61 Si exp (dVl) ' 

which, it is believed, results in a near-optimum solution sequence. 13 
The value of K varies from iteration to iteration and is different for 
every variable. 

Once the value of K t is determined, solution of (25) is accomplished 
as follows : 
Set 

D t = VI + KJ Si exp (dVl) + [1 + 6KJ Si exp (0^)]P,(7i) (27) 

where if, is picked according to equation (26). To obtain the new 
estimate, V) +l , set 

Yi = F* +1 (28) 

and establish an iterative procedure for solution of Y as 



yn+i _ Y" 



Y1 + K,I Bi exp (6Y n ) - D< 



(29) 



1 + eKJ Si exp (6Y •) 

where the superscript n indicates the iteration number in the scalar 
Newton-Raphson subloop. This is done for each element of [V]. 

A first estimate, F°, is formed as follows. 
Set 



then 



Z< = -±]n[K<I Si ]; 

if D, ;g Z, +0.08, F° = -D,, 

1 



and if Z), > Z { + 0.08, F- = - In D { + Z t . 

o 

The iterative procedure in Y f is terminated whenever 

I yn+i Y" 



m 



< e 



(30) 



(31) 



(32) 



where e is some suitably small constant, such as 10 -6 . 

The procedure outlined above, combined with the source stepping 
described in Section 3.3.2, provides an extremely powerful and rapid 
algorithm for solution of circuits with exponential nonlinearities. Solu- 
tion for most circuits is accomplished in very few iterations (fewer 
than 10), so that in combination with the efficient generation of the 
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equations, the overall method is competitive with analysis of linear 
systems. 

There is, however, one main limitation that must be dealt with as 
a special case. It arises from the formulation of the circuit equations 
and is discussed below. 



3.3.4 Treatment of Junction Cut-sets 

Whenever the circuit under consideration contains cut-sets of junc- 
tions, the matrix [A] of the system (1) becomes singular. This causes 
severe instability problems in the solution whenever the junctions in 
the cut-set are not (or barely) conducting. The problem, and its solu- 
tion, are best illustrated by an example. 

Consider the simple circuit in Fig. 1. The system of equations de- 
scribing this circuit is: 



1~] 






ll 




R 


v D1 


+ 


R 


[E]- 


1 






1 




fl_ 


V D2 




_ R- 





I 

'R 

"R 

for which the Jacobian is 



I sl [exp(6V m ) - 1] 
7 S2 [exp (dV D2 ) - 1] 



= 0, (33) 



[J] = 



~R ~ 6lsi exp ( 9Vd ^ 
R 



I 
~R 



-^ - dI S2 exp (6V D2 ) 



(34) 



which, by inspection, is seen to be extremely ill-conditioned whenever 
dig? exp (8V 1,2) ^ 1/2?. The conditioning of the system can be im- 
proved by subtracting the second equation from the first in (33) to 
vield: 



1^ 

R 



1 

R 


+ 






+ 







[E] 



7 sl [exp(0F m ) - 1] 



= 0, (35) 



_/ S2 [exp (6V D2 ) - 1] - 7 sl [exp(0V D1 ) - 1]_ 
for which the Jacobian is 



J = 



-| - 6I S1 exp (dV Dl ) 
dl si exp (6V m ) 



_1_ 
R 

-ei S2 exp (ev D2 )j 



(36) 
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I 



R 

-VW 






Fig. 1 — Junction cut-set example. 

The Jacobian shown in (36) can be scaled very simply to produce a 
well-behaved system. 

In the above example, the treatment of the cut-set problem was 
determined by inspection. Junction cut-sets can be found in any net- 
work by forming the so-called "L-tree" 16 of the network and forming 
the fundamental loop (or cut-set) matrix. The L-tree preference order 
dictates that junctions be made links with all other types of elements 
retaining their order. The treatment is simple, once junction/ current 
source cut-sets are found. For each cut-set, one row of the [4] and 
[B] matrices corresponding to one of the offending junctions is set 
to zero. The appropriate additions and subtractions in the vector 
[iV(V)] are performed, and the system scaled. A similar procedure is 
used by Shichman. 14 

3.4 Parameter Perturbation and Monte Carlo Analysis 

In addition to methods of solution of the nonlinear circuit equations, 
a method of statistically perturbing circuit parameters is required. 
The two are then combined in an overall strategy. 

3.4.1 Parameter Perturbation 

It is desired to generate correlated random variables with a fixed 
range, corresponding to the tolerance set by the user. Since large ar- 
rays of correlated numbers have to be generated, a parametric rep- 
resentation is used that correlates random variables by the use of 
"pivots." 1 In addition, a linear additive statistical model is used which 
ensures that parameters stay within tolerance. Generation of fixed 
interval correlated random numbers can be illustrated by a simple 
example: 

Assume x , x x and x 2 independent, each from a distribution of mean m 
and variance cj* . Two correlated random variables, y x and y 2 can be 
generated as 

ft = (1 — |Xi|)zi + \ix , .^ 

y 2 = (1 — |X 2 |)x 2 + \ 2 x , 
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where Ai, A a are tracking coefficients and x is serving as a "pivot." 
Note that if the x { are in the interval ( — 1, 1), then the y t will be in 
this interval also. The correlation factor p 12 is easily determined, 



P12 — 
resulting in 



E(y x y 2 ) - E(y,)E(y 2 ) ^ 



AiXa^io — XiX 2 Wlp 



(39) 



Pl2 " U(l " IXiDVi; + Xfaf.][(l - |X 2 |) 2 cr 2 , + Xjaf.]}* 
with 

E(3Ji) = (1 - |Xi|)m, + XiWo, 
E(y 2 ) = (1 — |X 2 |)ra 2 + X 2 m , 

2 /I l\ l\2 2 -v2 

0" yi = (1 — |Xi|) 0-,, — Ai<T Xo , 

o\, = (1 — |X 2 |) a Xl + X 2 o" Xo . 
For the special case of Xi = X 2 = X, <r x 2 o = a x \ — tr x \ = a and m = 0. 

X 2 
pi2 = 1 - 2 |X| + 2X 2 ' (40) 

The model actually used allows correlation to two pivots (or other 
parameters) as 

y { = (1 — |X,| - \i)i\)Xi + X.Xoi + r)iX 02 (41) 

and normalized random factors for the various parameters generated as 

U - 1 + t Vi (42) 

where t is the tolerance. 

The independent variables, x h are generated as follows: A piece- 
wise-linear probability density shape is supplied by the user in the 
form of a table in arbitrary units for distributions other than uni- 
form, normal or log normal which are "built in." This table is scaled 
and extended to include the cumulative density function which is a 
piecewise-quadratic on the interval ( — 1, 1). A random variable from 
a uniform distribution is generated by a Tausworthe random number 
generator 17 and transformed to the desired distribution by quadratic 
interpolation from the cumulative density function. Parameter per- 
turbations are then calculated according to equations (41) and (42). 
Nominal (design) values of parameters are taken to be the median 
value of their corresponding distributions. 
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3.4.2 Monte Carlo Analysis 

Parameter perturbations and analysis are combined in a standard 
tolerance analysis loop. The nominal solution is taken as the starting 
guess for each random design and the solutions along with the pa- 
rameter values stored on disk for later post-processing. Work is in 
progress to allow various other procedures (such as tuning or adjust- 
ment) in the loop. 

IV. IMPLEMENTATION 

The techniques and methods described in the preceding section al- 
low fairly efficient nonlinear statistical design. The details of their 
implementation can, however, spell the difference between success and 
failure in a general-purpose program, as well as the input/output 
features and structure of the program. Described below are some of 
the more critical aspects of the design and structure of a general- 
purpose nonlinear tolerance analysis program for IBM series 360 
computers. 

4.1 Memory Allocation 

As in any large program, there exists a conflict between efficient 
use of memory and speed of execution. In addition to the program 
itself, memory is required for the various coefficient matrices, the 
variables and outputs, as well as the various circuit parameters such 
as element values, model parameters, statistical data and topological 
information. 

The program itself, written largely in fortran-iv, is overlayed with 
major divisions separating (i) input and initial handling of data, 
(ii) generation of the various topological matrices, (Hi) generation of 
the equations, (iv) analysis, and (v) output. Communication among 
the various overlays is via a labeled common structure. Data at input 
time is handled largely via fixed-dimension arrays which are then 
partly overwritten by run-time data during execution. Run-time data 
is stored in a dynamic linear array with pointers used for addressing. 
This data includes such arrays as the A, B, C and D matrices, the 
variables, the Jacobian, various tabled quantities, and so on. In this 
way, efficient use is made of fast-access memory in that only data that 
is needed is stored. At the same time, this structure does not degrade 
the speed of execution. 

4.2 Algorithms 

The IBM 360 series is not well suited for applications such as described 
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here and special care is required in the implementation of the various 
algorithms in addition to standard good programming practices. All 
floating-point computation in the program is done in double precision 
(8 bytes). It was found that some of the matrix-handling subroutines 
had to be written in Assembly language in order to achieve any efficiency 
at all. For example, one routine that multiplies two matrices one of 
which has only 0, ±1 entries could be speeded up by a factor of 5-10 
by direct coding in Assembly language. Assembly language coding 
is also required for the equation solution and matrix inversion routines. 

Matrix sparsity is used to advantage in the solution of simultaneous 
equations as in equation (14) by column reordering 18 and stability 
is preserved by row-pivoting. Ill-conditioning is detected by moni- 
toring the magnitude of the smallest pivot used in the gaussian elim- 
ination process. 

Other small details in the programming are equally critical. It is 
of utmost importance in an application such as this to preserve as 
much numerical precision as possible since many mathematical steps 
are required before a solution is attained. For example, the analysis 
requires evaluation of quantities such as 7 s [exp (9V) - 1]. For values 
of V close to zero, a call to the exponential routine and subsequent 
subtraction of the constant 1 can yield an inaccurate result. For this 
situation, a series evaluation of the function is used. 

4.3 Data Reduction and Display 

Some care has to be taken in handling the voluminous data pro- 
duced by the program. Experience has shown that it is almost im- 
possible to determine beforehand how to analyze and display the 
output data, at least until a preliminary investigation of the results 
is available. In addition, any extensive Monte Carlo analysis of 
most circuits is likely to be expensive (despite the efficient algorithms) 
so that it becomes worthwhile developing a flexible post-processing 
scheme. For these reasons, the philosophy adopted here is the follow- 
ing: 

(i) All output data, including parameter values of every "im- 
portant" or expensive analysis, is stored in a permanent disk 
file for later access. This raw data may be later reduced or 
displayed in whatever way the user sees fit. 
(u) Extensive post-processing capability is available immediately 
following the analysis, accessed via the input language to the 
program. This facility allows scatter plots and histograms of 
any outputs or parameters to be produced on-line, including 
the printing out of extremal cases and various statistics. The 



1190 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1971 

facility also allows separation of the data by temperature or 
any other parameter. 
(in) Hard-copy printout of all data is produced whether or not the 
user requires it. In this way an expensive run is not lost even 
if the disk file is destroyed. Experience has also shown that 
sometimes it is desired to inspect the raw data weeks after 
initial analysis. The hard-copy output provides this facility. 

V. SAMPLE PROBLEMS 

Two sample problems are presented, both of which Monte Carlo 
analysis has verified to be good designs. The transistor model 2 used in 
both examples is shown in Fig. 2. In both problems silicon resistors 
(including the base and collector resistances of the transistor model) 
were allowed to vary ±15% within a normal distribution truncated 
at ±3o-, with resistor values tracking to ±5%. This is illustrated in 
the scatter diagram shown in Fig. 3. The intercept currents were picked 
from a log-normal distribution ranging from l/4/ So to 4/ So and cor- 
related by a factor of 0.85. The /3s of the devices were picked from a 
triangular distribution, typically ranging from l/2/3 to 20 o and cor- 
related by a factor of 0.3. 

The first example, a constant current source, used in a D/A converter, 
is shown in Fig. 4. Resistors Rl through R4 are thin film tantalum 
resistors with a tolerance of ±2%, with the exception of R4 which was 
assigned a tolerance of ±0.5%. R5 represents the load and was allowed 
to vary ±25%. Analysis of this circuit verified that the output current 
is essentially insensitive to all parameters except the power supplies 
and R4. A scatter plot of the output current versus the value of R4 
is shown in Fig. 5 for the case where the power supplies are held fixed. 
Figure 6 shows a histogram of the values of output current with all 
parameters varying. For this case, the twelve-volt supply was assigned 
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Fig. 2 — Transistor model. 
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Fig. 3 — Resistor correlation, range of each axis ±15%. 
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Fig. 4 — Constant current source. 
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Fig. 5 — Current source output (vertical axis) versus R4 (normalized) . 
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Fig. 6 — Current source output with all parameters varying. 
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to a tolerance of ±2%, with the six-volt supply tracking to ±1%. 
Analysis time for this example was approximately one second CPU 
time per random design on an IBM 360/65 computer. The great majority 
of this time represents overhead in the form of subroutine calls, various 
bookkeeping operations and writing data on disk and the printer. 
Solutions were carried to approximately seven digits of accuracy with 
each random design requiring typically 3-5 iterations. 

The second example shown in Fig. 7 is a silicon integrated opera- 
tional amplifier designed at Bell Laboratories. Of interest in this 
example is the output offset voltage with the inputs grounded. In 
addition, minimum and maximum (worst case) currents in the col- 
lectors of T13, T14 and T16, as well as dc gain were sought. This 
circuit represents a rather large simulation with 48 junctions, 76 resis- 
tors, 79 nodes and 127 branches. Figure 8 shows a histogram of the 
output voltage at room temperature with the power supplies held fixed. 
The range of offsets for this circuit was found to be slightly higher 
than predicted by approximate hand analysis. Analysis time for this 




Fig. 7 — Operational amplifier. 
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Fig. 8 — Operational amplifier output offset voltage. 

example was 8 seconds CPU time per statistical design, some of this 
time being bookkeeping overhead. It is expected that this time will 
be cut down considerably with some reprogramming. As in the pre- 
vious example, typically 3-5 iterations were required per statistical 
design for a seven-digit accuracy in junction voltages. 

VI. CONCLUSIONS 

Analysis techniques and programming considerations have been 
presented that allow reasonably economical tolerance analysis of 
nonlinear "dc" circuits of reasonable size. Many of the ideas presented 
have evolved from past experience with ac tolerance analysis and will 
most probably be modified as experience with nonlinear statistical 
design becomes more plentiful. It is already apparent, however, that 
the trend in the near future will be to larger scale integration of 
circuits for which some of the present analysis techniques will likely 
be inadequate. Research is in progress in analysis methods capable 
of coping with large and complex circuits, as well as methods to make 
the present techniques even more efficient. 
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